
suppressPackageStartupMessages({
  library(sf)
  library(terra)
  library(exactextractr)
  library(dplyr)
  library(readr)
})

# ---------- Paths ----------
if (!requireNamespace("here", quietly = TRUE)) install.packages("here", quiet = TRUE)
base_dir <- normalizePath(here::here(), winslash = "/")

SUBPLACE_FILE <- file.path(base_dir, "Data", "Shape_Files", "SP_SA_2011", "SP_SA_2011.shp")
subplace_file <- SUBPLACE_FILE
raster_file   <- file.path(base_dir, "Data", "Remote_Sensed", "GHS_POP_100_2020.tif")

out_dir <- file.path(base_dir, "Processed_Data", "Remote_Sensed")
dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)

out_csv <- file.path(out_dir, "subplace_population_ghsl100m.csv")
out_rds <- file.path(out_dir, "subplace_population_ghsl100m.rds")

stopifnot(file.exists(subplace_file), file.exists(raster_file))

# ---------- Read ----------
sp <- st_read(SUBPLACE_FILE, quiet = TRUE)
# Fix missing CRS: assume WGS84 if undefined
if (is.na(st_crs(sp))) st_crs(sp) <- st_crs(4326)
sp$SP_CODE <- as.character(sp$SP_CODE)

# Read raster  (fix: use raster_file)
r <- terra::rast(raster_file)
r_crs <- terra::crs(r, proj = TRUE)  # WKT/proj string
cat("Raster CRS:", r_crs, "\n")

# ---------- Reproject subplaces to raster CRS for extraction ----------
sp_r <- st_transform(sp, r_crs)
sp_r <- st_make_valid(sp_r)
sp_r <- sp_r[!st_is_empty(sp_r), , drop = FALSE]

# ---------- Area in an equal-area meter CRS ----------
# If the raster CRS isn't projected in meters, use world Mollweide for area only.
is_longlat <- sf::st_is_longlat(sp_r)
area_crs <- if (is_longlat) {
  sf::st_crs("+proj=moll +lon_0=0 +datum=WGS84 +units=m +no_defs")
} else {
  sf::st_crs(sp_r)
}
sp_area <- st_transform(sp_r, area_crs)
area_m2  <- as.numeric(st_area(sp_area))
area_km2 <- area_m2 / 1e6

# ---------- exactextractr (coverage-weighted) ----------
extract_df <- exact_extract(
  r,
  sp_r,
  fun = function(values, coverage_fractions) {
    keep <- !is.na(values) & coverage_fractions > 0
    if (!any(keep)) return(data.frame(pop_density_100m_mean = NA_real_, population = 0))
    v <- values[keep]; w <- coverage_fractions[keep]
    data.frame(
      pop_density_100m_mean = sum(v * w) / sum(w),  # mean people per 100m pixel
      population            = sum(v * w)            # total people
    )
  },
  progress = TRUE
) |>
  bind_rows()

# ---------- Assemble & write ----------
out <- sp_r |>
  st_drop_geometry() |>
  transmute(
    SP_CODE,
    area_sp_m2  = area_m2,
    area_sp_km2 = area_km2
  ) |>
  bind_cols(extract_df) |>
  mutate(population = ifelse(is.na(population), 0, population)) |>
  select(SP_CODE, pop_density_100m_mean, area_sp_m2, area_sp_km2, population)

# Write outputs (fix: use out_csv / out_rds)
write_csv(out, out_csv)
saveRDS(out, out_rds)

